#include "mcRandVector.h"

// construction & destruction
irmm_RandVector::irmm_RandVector()
{
	rf = NULL;
}

irmm_RandVector::~irmm_RandVector()
{
	if (rf!=NULL)
		delete []rf;
	rf = NULL;
}

void irmm_RandVector::Set(mcmcVector<double> &mean, mcmcMatrix<double> &sigma)
{
	int		i;

	// copy mean and covariance matrix
	mx_mean = mean;
	mx_sigma = sigma;

	// compute eigen values and eigen vectors
	mx_sigma.Eigen(mx_eigen_value,mx_eigen_vector);

	if (rf!=NULL)
		delete []rf;
	rf = new RF_uniform[mx_eigen_value.rows()];
	for (i=0; i<mx_eigen_value.rows(); i++)
		rf[i].Set(-sqrt(fabs(mx_eigen_value(i))),sqrt(fabs(mx_eigen_value(i))),1.0);
}

mcmcVector<double> irmm_RandVector::Sampling(void)
{
	mcmcVector<double>	mx_samples;
	int				i;

	mx_samples.SetRowVector(mx_eigen_value.rows());
	for (i=0; i<mx_eigen_value.rows(); i++)
		mx_samples(i) = rf[i].Sampling();
	mx_samples = mx_mean + mx_eigen_vector * mx_samples;

	return mx_samples;
}

double irmm_RandVector::GetPdf(mcmcVector<double> &mx_value)
{
	mcmcVector<double> mx_temp;
	double			dpdf;
	int				i;

	mx_temp = mx_value-mx_mean;
	mx_temp = mx_eigen_vector.T()*mx_temp;
	dpdf = 1.0;
	for (i=0; i<mx_value.rows(); i++)
		dpdf = dpdf*rf[i].pdf(mx_temp(i));
	return dpdf;
}

void irmm_RandVector::Diffuse(mcmcVector<double> &mx_value)
{
	mcmcVector<double>	mx_temp;
	double			temp;
	int				i;

	// transfer in the eigen space
	mx_temp = mx_eigen_vector.T()*(mx_value-mx_mean);
	temp = 0.0;
	for (i=0; i<mx_temp.rows(); i++)
		temp += mx_temp(i)*mx_temp(i)/mx_eigen_value(i);
	if (temp>1.0)	// outside the ellipsoid
	{
		// find the intersection with the ellipsoid
		temp = sqrt(1.0/temp);
		for (i=0; i<mx_temp.rows(); i++)
			mx_temp(i) = mx_temp(i)*temp;
		mx_value = mx_mean+mx_eigen_vector*mx_temp;
	}
}
